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Abstract — Orthogonal Frequency Division Multiplexing 
(OFDM) is the key component of many emerging broadband 
wireless access standards. The resource allocation in OFDM 
uplink, however, is challenging due to heterogeneity of users' 
Quality of Service requirements, channel conditions, and 
individual resource constraints. We formulate the resource 
allocation problem as a non-strictly convex optimization 
problem, which typically has multiple global optimal solutions. 
We then propose a reduced primal-dual algorithm, which is 
distributed, low in computational complexity, and probably 
globally convergent to a global optimal solution. The performance 
of the algorithm is studied through a realistic OFDM simulator. 
Compared with the previously proposed centralized optimal 
algorithm, our algorithm not only significantly reduces the 
message overhead but also requires less iterations to converge. 

I. Introduction 

Orthogonal Frequency Division Multiplexing (OFDM) is a 
promising technology for future broadband wireless networks. 
In OFDM, the entire frequency band is divided into a large 
number of subchannels, and network resource can be allocated 
flexibly over each of the subchannels. Because of this, OFDM 
enjoys several key advantages such as robustness against 
intersymbol interference and multipath fading, and no need 
for complex equalizations. 

In this paper, we consider the resource allocation in a 
single cell OFDM uplink system, where multiple end users 
transmit data to the same base station. This is motivated by 
several practical wireless systems, such as WiMAX/802.16e, 
LTE for 3GPP, and UMB for 3GPP2. Given the channel 
conditions of the users at a particular time, we need to 
determine which subset of users to schedule (i.e., allowed to 
transmit with positive rates), how to match the subchannels 
with the scheduled users, and the power allocation across these 
subchannels. 

We formulate the resource allocation problem as a weighted 
rate maximization problem, which is motivated by the 
gradient-based scheduling framework in Q~|-|[3). This problem, 
however, is quite challenging to solve due to the heterogeneity 
of users' Quality of Service requirements, channel conditions, 
and individual resource constraints. In this paper, we propose 
the first distributed primal-dual algorithm in literature that 
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achieves the optimal resource allocation in uplink OFDM 

systems. Our key contributions are: 

. Realistic OFDM model: we consider several important 
practical constraints such as self-noise and per subchannel 
per user SNR limits. These constraints are typically 
ignored in previous work but have significant impacts on 
the optimal solution as well as algorithm design. 

• Optimal algorithm with global convergence: the proposed 
algorithm is provably globally convergent to one of 
the global optimal solutions of the resource allocation 
problem, despite of non-strict convexity of the problem 
under which setting primal-dual algorithms may not be 
able to converge (e.g., Pfl-ll6l). 

. Distributed algorithm with scalable performance: the 
proposed algorithm is distributed, requires simple local 
updates, demands only limited message passing, and is 
highly scalable with network size. 

• Simple algorithm with fast convergence: the proposed 
algorithm only updates a subset of all decision variables. 
Compared with the previously proposed centralized algo- 
rithm, our algorithm reduces complexity and requires less 
iterations to converge. 

In Section [XT] we present the network model and problem 
statement. The reduced primal-dual algorithm is proposed 
in Section [HI] and its convergence behavior is analyzed in 
Section |IV] Simulations results on convergence and optimality 
are given in Section [V] and we conclude in Section [VT] 

A. Related Work 

Most previous work on resource allocation in OFDM sys- 
tems focused on the downlink case, where the base station 
sends traffic to multiple end users with a total power constraint. 
Compared with the uplink case, the optimization problem in 
the downlink case is easier to solve and a centralized algorithm 
is reasonable to implement (e.g., Q). Due to different resource 
constraints, however, the algorithms proposed for the downlink 
case can not be directly applied to the uplink case. 

Uplink OFDM resource allocation only receives attention 
recently, e.g., lf8l- lfT4l . In (8), the problem was formulated in 
the framework of Nash Bargaining with a focus of fair resource 
allocation. The authors of |9] proposed a heuristic algorithm 
that tries to minimize each user's transmission power while 
satisfying the individual rate constraints. In iflOl . the author 
considered the sum-rate maximization problem and derived 
algorithms based on Rayleigh fading on each subchannel. 
The authors in |[TTI - lfT4l proposed several heuristic algorithms 
to solve a problem similar as the one considered here with 
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additional integer channel allocation constraints. None of the 
previous literature focused on solving the uplink resource 
allocation problem optimally. 

In a recent work ||T5l . the authors proposed a dual -based 
centralized algorithm to solve the resource allocation problem. 
In the centralized algorithm, the base station needs to col- 
lect the complete network information (e.g., users' weights, 
resource constraints, and uplink channel conditions) before 
the algorithm starts. It also needs to have sophisticated com- 
putational capability to perform the dual-based algorithm. 
Finally, recovering the primal optimal solution after the dual 
algorithm converges is not trivial. Our algorithm overcomes all 
three bottlenecks since it is distributed, simple, and provably 
converges to the optimal primal variables. 

This paper lies in the area of using primal-dual algorithms 
for non-strictly convex optimization problems. Under the non- 
strict convex setting, primal-dual algorithms in general may 
not converge to the optimal solution, as observed in litera- 
ture j4], J6). In j6), the authors characterized the convergence 
behavior of a certain primal-dual algorithm in a P2P setting, 
and gave a sufficient condition under which the algorithm 
is guaranteed to converge. Our analysis is motivated by the 
convergence proof in J6), but we deal with a different set of 
challenges in a totally different problem setting. 

This work is a generalization of our previous work 0, 
where we presented some preliminary results on the primal- 
dual algorithm for a simple OFDM model. Here we consider 
a much complicated model including various constraints such 
as self-noise and SNR limitations. Moveover, we design a 
reduced primal-dual algorithm by considering the relationship 
that the variables should satisfy at the optimal solution. In 
other words, the dynamics of our newly proposed algorithm 
is constrained by a high dimension nonlinear manifold. All 
the above require different techniques in terms of proving 
the optimality and convergence. We also perform extensive 
simulations to illustrate the performance of the proposed 
algorithm together with comparison with the previous dual- 
based centralized algorithm. 

II. Problem Statement 

We consider a single OFDM cell, where there is a set 
M = {1, ... , M] of users transmitting to the same base station. 
Each user i e M has a priority weight w/QThe total frequency 
band is divided into a set N = {!,... ,N] of subchannels 
(e.g., tones/carriers). A user i e M can transmit over a subset 
of the subchannels, with transmit power py over subchannel 
j e N and the total power upper bound, i.e., £y Pij ^ Pi- F° r 
channel j, it is allocated to user i with fraction jty, and the 
total allocation across all users should be no larger than 1, i.e., 

We define gy as the received signal-to-noise ratio (SNR) per 
unit power for user i on subchannel j. As in Q, this model 
can also incorporate various sub-subchannelization schemes 
where the resource allocation is performed in terms of sets of 

'We will later discuss how these weights are derived in practice. 



TABLE I 
Key Notations 



Notation 


Physical Meaning 


AJ 
1\ 


total number of subchannels 


\f 
J\ 


set of all subchannels 


M 


total number of users 


\A 
JV\ 


set of all users 




user index 


i 


subchannel index 


Wi 


user Z's (dynamic) weight 


eij 


normalized SNR on subchannel j for user i 


PU 


power allocated on subchannel j for user i 


Xij 


fraction of subchannel j allocated to user i 


Pi 


maximum transmit power for user i 


Si) 


maximum SNR coefficient on subchannel j for user i 


P 


self-noise coefficient 



frequency bands in the frequency domain or with a granularity 
of multiple symbols in the time domain. We further assume 
that the channel conditions do not change within the time of 
interests, i.e., we are looking at a resource allocation period 
smaller than the channel coherence tirn^l 

The key notations used throughout this paper are listed in 
Table U We use bold symbols to denote vectors and matrices 
of these quantities, e.g., w = {w,,V/}, e = {gy, Vz,j}, p = 
[Pij, W, j], and x = {x u , V/, ;'}. 

Our objective is to maximize the weighted sum of the users' 
rates over the feasible rate-region defined as follows, 

*•>-{-««?:■>- §W 1+ 3^l v,e H 

(i) 

where r, is the achievable rate of user (J and (x,p) e X are 
chosen subject to 

J]x u <l,VjeN, (2) 

i 

J^Pij < P b Vi 6 M, (3) 

j 

and the set 

X := l(x,p) > : < x u < 1,0 < py < —. //../]. (4) 
I ey I 

Here, fl represents the level of "self-noise" because of imper- 
fect carrier synchronization or inaccurate channel estimation 
(e.g. HU, IfTTl ). sy is a maximum SNR constraint on sub- 
channel j for user i, modeling for example limited choices of 
available modulation and coding schemefl 

In many OFDM standards, xy is constrained to be an inte- 
ger, in which case we add the additional constraint xy e {0, 1 } 
for all j. The integer constraint makes the resource allocation 
very difficult to solve, and various heuristic algorithms dealing 
with such constraint are proposed in |[T2l - lfT5l . In this paper, 

2 This is particularly suitable for fixed broadband wireless access (part of 
the IEEE 802.16 standard), where users are relatively static. 

3 For notation simplicity we normalize the bandwidth to be 1 during the 
analysis. A realistic value will be considered in the simulations (Section IvY 

4 We restrict our attention to the interesting "resource constrained" case, 
where constraints {5) are tight for all users at an optimal solution with the 
maximum SNR constraints taken into consideration. 
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we will ignore this integer constraint and focus on the rate 
region defined by © to ©. The corresponding solution 
typically contains fractional values of xy's. There are several 
practical methods of achieving these fraction allocations. For 
example, if resource allocation is done in blocks of OFDM 
symbols, then fractional values of Xy can be implemented by 
time-sharing the symbols in a block. Likewise, if the number 
of subchannels are large enough so that the channel conditions 
do not change dramatically among adjacent subchannels, then 
the fractional values of x/fs can also be implemented by 
frequency sharing (e.g., IfTSl ). 

In this paper, we focus on solving the following problem 



max > Wirt, 

reR(e) ' 



(5) 



where rate r, and rate region 7?(e) are given by ©. 

Finally, we note that the priority weights w,'s are motivated 
by the gradient-based scheduling framework presented in UJ- 
© . Let's assume each user ; has a utility function t//(W/,,) 
depending on its average throughput W,,, up to time t. It has 
been shown that in order to maximize the total network utility 
Yui Uj(Wjj) in the long run, it is enough to solve Problem Q 
during each time slot t with w-, = dU;(W- u )ldW;j. 

III. A Reduced Primal-Dual Algorithm 

We can rewrite problem © in variables x and p as follows: 
Problem 1 (Weighted Rate Maximization): 



max 



ieM 



log 1 + 



Xij +Ppue k 



(6) 



subject to the per subchannel assignment constraints in <J2J and 
the per user power constraints in ©. Set X is given in (0J. 

Although the objective function in © is concave, its deriva- 
tive is not well defined at the origin (x — Q,p — 0). This 
motivates us to look at the following e-relaxed version of 
Problem [TJ 

Problem 2 ( e-relaxed Weighted Rate Maximization ): 



max 

(x.p)eX 



^ Wi ^(Xij + etj) log II + - 

ieM jeN \ X ' 



PU e U 



i+fiPijeij + Eij 



(7) 



where constants ey take small positive value for all i and j. 
The constraint set remains the same as in Problem [TJ 

By such relaxation, the objective function in © now has 
derivative defined everywhere in the constraint set X. Thanks 
to the continuity of the objective function, the optimal value 
to Problem [2] can be arbitrarily close to that of Problem [TJ if 
e = [eij, Vi, j] is chosen to be small enough. 

The constraint set of Problem [2] is convex, and the objective 
function in © is continuous and non-strictly concave. As such, 
Problem [2] has multiple optimal solutions, and there is no 
duality gap between it and its dual problem. 

The existence of derivatives allows us to write down a 
primal-dual algorithm to pursue the optimal solution to Prob- 



lem |2] The Lagrangian for Problem [2] is as follows, 



L(A, fi, x,p) := V Wiixtj + eij) log(l + 1^ ) 

Xj + Pptjeij + e t j 

+ 2 Ai(Pi - 2 PiJ ) + 2^(1-2 *y)- (8) 

< j j 

By strong duality theorem, the optimal primal and dual solu- 
tions must satisfy KKT conditions, i.e., for all i and j, 

tij>0, J]xij<l, x t j -1) = 0, (9) 

i i 
At > 0, Yj PH P i> PiJ - P i) = °. ( 10 > 

Xij Si j 



J J 

Xij > 0, < pij < 



(11) 



Xij(Jij{Xij, pij)-pj)<0, (12) 
Pij(gij( X U> Pij) — ^i)< 0, (13) 

where fij(-) and gij(-) are gradients of the objective function 
in © with respect to x t j and p t j, respectively, and are given 
by 



fiMipPij) = Vf;l0g(l + 



PU e U ) 
xtj +/3pijeij + e tj 
w^Xij + e^pijetj 



(x u +/3pijeij + ejj)[xij + (J3 + 1 )/>;/•,, + e,y] ' 



and 

gij(Xij,Pij) 



wtetjixij + etj) 2 



(x^ +Ppijeij + eij)[xjj + (J3+ \)pnc,; + <7.<l 



The last two KKT conditions in (fT2l and (fTTt become 
equalities if x, y > and p t j > 0, respectively. It can be verified 
that the optimal solutions of Problem [2] satisfying above KKT 
conditions, are exactly the saddle points of the Lagrangian 
function in ©. Since the primal problem has at least one 
solution, the saddle point exists. 

For notation simplicity, we define (a) + = max(a, 0), 



and 



a, b > 0, 

max(a, 0), otherwise, 



< b < c, 



(a)$ c = { max(a,0), b < 0, 



min(a, 0), otherwise. 



To pursue saddle points of the Lagrangian function, we start 
by considering a standard-form primal-dual algorithm, then we 
derive a new reduced primal-dual algorithm which has less 
complexity and faster convergence speed than the standard 
one. 
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Motivated by the work in Q, we consider the following 
standard-form primal-dual algorithm: V/, j, 



Xij 


= ^(./n'V;/. /;;/)-//;); . 


(14) 


Pij 


= k ij(gij( x U> Pij) - A ifp l} ;JBL ' 


(15) 


l-'j 


i 


(16) 


Ai 


= kf(J]pij-Pi)X> 


(17) 



where kf^k^kf^ and kj are constants representing adaption 
rates. Here the derivatives are defined with respect to time. 

We call a point (x,p, fx, A) an equilibrium of the standard- 
form algorithm if and only if the corresponding derivatives in 
(fT4l to ( fTTI ) are zero for all i and j. We can show that the 
set of equilibria of the standard-form primal-dual algorithm is 
equivalent to the set of global optimal solutions of Problem [2] 

In the standard-form primal-dual algorithm, all variables x^, 
Pij, fij, and A/ are dynamically adapted. This might lead to 
performance concerns in terms of high complexity and slow 
convergence. One way to address this concern is to reduce 
the number of dynamically adapting variables. We achieve this 
goal by constraining the algorithm trajectories onto a manifold 
which includes all optimal primal and dual solutions. 

We study the following manifold by setting < fT3T > to zero, 
i.e. ViJ, 



= (Hij(xi h p l ,)-A,y „j, tJ , 

Pij' jr 



which in turn implies 



gij(xij, P ij) = A h if <Pij <*BL, 
gij(xij,pij) < Ai, if p u = 0. 



(18) 



(19) 



Clearly, the optimal primal and dual solutions must lie on the 
above manifold. 

After simplification we get the following expression of the 
manifold: 

XjjSij 



Pij = mm 



en 



(20) 



where h{j is denoted by 

^/l+W+l)^-(2 J 6+l) 



2(303 + l)e, 



(Xij + €ij)) 

'WiCij-Ai 



when {3 + 0, and 



(Wtejj-Ai \+ 
hij= — (Xij + eij)) 

V A-fiij >K:,:-. 



when f3 — 0. 

Substituting (1201 into ( TT4T > to (fTTb . we obtain a new reduced 
primal-dual algorithm as follows: 



Algorithm RPD: Reduced Primal-Dual Algorithm 

X U = %j(fij( x ij>Pij)-P-jfx tj > 



Pij = 



Pu - Pi) 



+ 

.1, • 



kj 



(21) 
(22) 

(23) 
(24) 



Proposition 1: The set of equilibria of Algorithm RPD is 
the same as the set of global optimal solutions of Problem [2] 

This means that if the reduced primal-dual algorithm con- 
verges, it reaches a global optimal solution of the e-relaxed 
weighted rate maximization problem. 

Compared to the standard-form algorithm in (fT4l to ( fTTI ) in 
which p is dynamically adapted, p in the new Algorithm RPD 
is directly computed from x and A. Consequently, the reduced 
Algorithm RPD has less dynamically adapting variables than 
the standard-form algorithm, and hence is less complex and is 
expected to converge faster. 

Similar as the standard-form algorithm, the reduced Algo- 
rithm RPD can also be implemented in a distributed fashion 
by end users and the base station. A end user i is responsible 
of updating x,/s as well as dual variable Aj locally. During 
each iteration, it sends the latest values of x,/s to the base 
station, but not the py's or Aj. The base station is responsi- 
ble for updating dual variables fij's for all subchannels and 
broadcasting to the users. In particular, the base station does 
not need to know users' channel conditions, power constraints, 
or priority weights. Both the communication complexity and 
computation complexity per iteration are 0(MN). 

Despite of its advantage, we need to first show that trajec- 
tories of Algorithm RPD converge to its equilibria before we 
can apply it to solve Problem [2] 

IV. Convergence of the Reduced Primal-Dual Algorithm 

In Algorithm RPD, we remove power allocation variables 
Pij's from the dynamically adapting variables. Now, the dif- 
ficulty is how to ensure the global convergence of Algorithm 
RPD. The key challenge of the convergence proof is the non- 
strict concavity of the objective function in (|7). It has been well 
observed in the literature that although primal-dual algorithms 
can globally converge to the optimal solution of a strictly 
concave optimization problem, they may oscillate indefinitely 
and fail to converge when applying to a non-strictly concave 
optimization problem, even setting small enough update step- 
sizes El-El, fl9l . 

In this section, we successfully prove the sufficient condi- 
tions under which Algorithm RPD can globally and asymptot- 
ically converge to a global optimal solution of Problem [2] The 
proof is organized as follows. We first show in Theorem Q] that 
the trajectories of Algorithm RPD converge to an invariant set 
containing all global optimal solutions of Problem [2] Then in 
Lemma Q] we show that the trajectories of Algorithm RPD 
might form limit cycles within the invariant set if this set 
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contains non-optimal solutions. In Theorem [2] we present the 
conditions under which the invariant set only contains global 
optimal solutions of Problem [2] without limit cycles. Finally 
we show in Corollary Q] how to choose the update stepsizes to 
satisfy the sufficient conditions in Theorem [2] 

Theorem 1: All trajectories of Algorithm RPD converge to 
an invariant set Vo globally and asymptotically. Furthermore, 
let (x* , p* , p* , X*) be a global optimal solution of Problem [2] 
and {x,p, p,\) be any point in set Vo, the following is true 
for all i and j, 

1) (x*,p*, p*, A*) is contained in Vq\ 

2) fij is nonzero only if 2***; = 1; 

3) Yij P*j — Pi' an d i s a positive constant; 

4) over set V , /;,( v ;; ./),,) = f,,{ \;,. />;,) = [I*, and 
gij(xij,pij) = gij(x*j,p*j) = A*; 



5) 



The proof of Theorem Q] can be found in Appendix [A] 
Results [TJ to |4]i will be used in later analysis. 

Result|5]l of Theorem[T]is of independent interest. It implies 
that although there can be multiple global optimal solutions to 
Problem|2j the effective SNR achieved by user i on subchannel 
j is the same in all solutions. 

Although all trajectories of Algorithm RPD may converge 
to the desired equilibria in Vo, they may also converge to non- 
equilibrium points in Vo (if there are any). We now study the 
conditions for Vo to contain only the desired equilibria, under 
which Theorem \T\ guarantees the convergence of Algorithm 
RPD to a global optimal solution of Problem [2] 

Plugging result |4|i of Theorem Q] into Algorithm RPD, and 
recalling that M is the total number of users and N is the total 
number of subchannels, we find that Vo is exactly the set that 
contains all trajectories of the following linear system in (|25| > 
to d27j over set {x > 0, p > 0). 

x = K x A\p* -K x A T x p, (25) 
K^Aix-Kn, (26) 
K A A 2 B{x + e) - K A P = 0, (27) 



A* 
A 



where K x is an MNxMN diagonal matrix with diagonal terms 
equal to fey's, K M is an N x N diagonal matrix with diagonal 
terms equal to A^'s, and K A is anMxM diagonal matrix with 
diagonal terms equal to kf's. B is an MNxMN diagonal matrix 



given by B — diag(bij, Vi, 7), where b-,j — ,'' . The matrix 

Ai = [In, ■ ■ ■ ,In] and has a dimension of N x MN, where 1^ 
is an identity matrix with dimension N. The matrix An has a 
dimension of M x MN and is given by 



llxAT 
















where li x ^ is an all one vector with dimension 1 by Af. 

For the above linear system, we have the following obser- 
vations. 



Lemma 1: For the linear system in d25l l to 1271 . we have 

1) every order Lie derivative of A^Bx is constant, that is 
Vn > 0: 

d" 

— AiBx = constant, 
dt" ' 

where t denotes time; 

2) starting from a non-equilibrium point, trajectories of 
x and fi, following d25T l and (|26| | respectively, do not 
converge and form limit cycles. 

The proof of Lemma Q] can be found in Appendix [B] 

Result [2] in Lemma Q] indicates that if (x, fi) starts from a 
non-equilibrium point, they will keep oscillating and never 
converge. Therefore, one way to guarantee the system conver- 
gence is to make sure that the invariant set Vo only contains 
equilibria points (i.e., global optimal solutions). Result Q] in 
Lemma Q] states that every order Lie derivative of AiBx is 
constant. By linear system theory, if the system state fj, is 
completely observable from A2BX, then constant A2BX will 
lead to (1 equal to and being constant. When is constant, 
x is constant according to (l25l l. Combining with the constraint 
that x > and A^Bx is constant, we can show that x is also 
zero if (1 is zero over the set {x > 0, p, > 0). 

In the following theorem, we state conditions for p to 
be completely observable from A2BX, and summarize its 
consequence on convergence of Algorithm RPD. 

Theorem 2: All trajectories of Algorithm RPD converge 
globally and asymptotically to the system equilibria if the 
following condition holds: 



A 2 BK X A\ 
K^AiK^Aj — crl 



has rank Af, 



(28) 



where <x denotes any eigenvalue of matrix K t 'A\K x A\. 

The proof of Theorem [2] can be found in Appendix [C] 

For the problem we studied in this paper, we can choose 
properly the update stepsizes of the algorithm to satisfy the 
condition in ( f28T >. 

Corollary 1: The condition (f28b in Theorem [2] is satisfied 
if both of the following are true 

. K x = kl (diagonal terms of K x take the same value k); 

. all diagonal elements of K M take different values. 
The proof of Corollary Q] can be found in Appendix [D] 

In this section, we have investigated the convergence of 
Algorithm RPD by combining both La Salle principle from 
nonlinear stability theory and complete observability from 
linear system theory. The proof shows that Algorithm RPD 
can globally and asymptotically converge to one of the global 
optimal solutions of Problem [2] when satisfying the conditions 
in Corollary Q] The convergence proof of Algorithm SPD is 
similar and will not be presented here due to space limitation. 

V. Simulation Results 

We show the convergence and optimality of the reduced 
primal-dual algorithm in Algorithm RPD over a realistic 
OFDM uplink simulator. Each user's subchannel gains e,y's 
are the product of two terms: a constant location-based term 
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Channel conditions of four users 



Equal weights 




5 10 
Channel index (0 to 15) 

Fig. 1. Channel conditions of 4 users and 16 channels 

picked using an empirically obtained distribution, and a fast 
fading term generated using a block-fading model and a 
standard mobile delay-spread model with a delay spread of 
10jusec. The system bandwidth is 5MHz consisting of 512 
OFDM tones, which is further grouped into 64 subchannels[| 
The symbol duration is 100yusec with a cyclic prefix of 10yusec. 

Unless otherwise specified, we assume the following pa- 
rameter setting throughout all simulations. The variables are 
initialized as = 1/M, py = P//N, fij — 0, and Aj 
0.01 x max/vt^ey) for all i and j. The update stepsizes in 
Algorithm RPD are chosen as £0 = 10~ 2 , = 10 _1 + Sj, and 
kf - 10~ 2 for all i and j. Here, s/s for all channels are chosen 
to be very small values and diverse from each other, in order to 
meet the requirement of Corollary Q] as the sufficient condition 
for system convergence. Each user has a total transmission 
power constraint P, = 2Watts. Users' channel conditions are 
randomly generated from the simulator, and users have equal 
weights Wi = 1 for all i, 

A. Optimal Resource Allocation 

We first demonstrate how the power and subchannel are 
allocated once Algorithm RPD converges to the optimal so- 
lution. In order to show the results effectively, we consider 
a small network with 4 users and 16 subchannels. Larger 
networks with be simulated in Sections [V-B| to |V-Dl The initial 
values of Aj's are chosen to be 0.1 x max ; -(w,-ey) for all ;. The 
self-noise coefficient is /3 = 0.01, and there are no maximum 
SNR constraints (i.e., sy = oo for all i and j). Fig. Q] shows 
the channel conditions of 4 users. The horizontal axis is the 
channel index (from to 15), and the vertical axis is the 
normalized channel gain <?y. It is clear that user 1 has the 
best channel condition and user 2 has the worst. We simulate 
Algorithm RPD under two different priority weight settings: 
(a) equal weights where w,- = 1 for each user i, and (b) unequal 
weights where w\ — 2, W2 — W3 — 1, and w\ — 0.5. 

5 Every 8 adjacent tones are grouped into one subchannel. This corresponds 
to the "Band AMC mode" of 802.16 d/e and can help to reduce the feedback 
overhead. For discussions on various ways of subchannelization, see (7J. 




Fig. 



Channel index (0 to 15) 
2. Optimal channel allocation under Algorithm RPD 



Equal weights 




Channel index (0 to 15) 
Fig. 3. Optimal power allocation under Algorithm RPD 



Equal weights 




Channel index (0 to 15) 
Fig. 4. Optimal rate allocation under Algorithm RPD 

Figures [2] [3] and [4] show the optimal channel allocation, 
power allocation, and rate allocation, respectively. In all three 
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Fig. 5. Case 1 for 40 users: primal and dual variable convergence of 
Algorithm RPD 



400 
£• 200 



Total weighted rate 



- Dual value (upper bound) 

- Primal feasible value (lower bound) 



son 



1000 
Relative errors 



1500 



200(1 



g 10 

> 

J 10 
I 




son 



1000 

Number of iterations 



1500 



2000 



Fig. 6. Case 1 for 40 users: total weighted rate convergence of Algorithm 
RPD 

figures, the top and bottom subfigures correspond to equal 
weights and unequal weights, respectively. Now consider a 
user with a good channel condition and a higher weight, e.g., 
user 1 . It is clear that user 1 is allocated with larger fractions 
of channels (Fig. |2j. Since the total power constraint for all 
users is the same, user 1 actually transmits with less power on 
average over the subchannels allocated to him (Fig. [3}. The 
net result is that he achieves a better data rate (Fig. |4). 

B. Algorithm Convergence 

1) 40 Users: Next we show the convergence of Algorithm 
RPD with 40 users and 64 subchannels. Here we assume 4 
cases as shown in Table HH 

For Case 1, Figure [5] shows the convergence of dual 
variables (upper two subgraphs, A, for 40 users and p.j for 
64 subchannels) and primal variables (lower two subgraphs, 
Yjj Pij for 40 users and Yii x U f° r 64 subchannels). 

In Fig. |6j the upper subgraph shows how the dual value 
and primal feasible value change with iterations. The dual 
value is an upper bound of the global optimal solution of 



TABLE II 

Four Cases 



Case 


1 


2 


3 


4 


Self-noise coefficient /? 








0.01 


0.01 


SNR constraints sn (dB) 


CO 


20 


oo 


20 



Dual variables I 



Dual variables n 




500 1000 1500 2000 
Total power allocation of all users 





500 1000 1500 2000 



Total channel allocation of all channels 
2 ii 



500 1000 1500 2000 
Number of iterations 



1.5 
1 

0.5 




500 1000 1500 2000 
Number of iterations 



Fig. 7. Case 2 for 40 users: primal and dual variable convergence of 
Algorithm RPD 
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Case 2 for 40 users: total weighted rate convergence of Algorithm 



Problem [2] The primal feasible value is a lower bound and 
is calculated as follows: given the primal values of p(t) and 
x(t) at iteration t, normalize so that they are feasible and the 
resources are fully utilized (i.e., pij(f) = Pij(t)Pi/(Tjj Pij(t)) 
and Xjj(t) = Xij(t)/{Yji Xij(t)) ), and calculate the achievable rate 
accordingly. The bottom subfigure shows the relative errors of 
two curves plotted in the upper subfigure. If we define the 
stopping criterion to be the relative error less than 5 x 10 -3 , 
then Algorithm RPD converges in 364 iterations. 

We also simulate Cases 2-4 in Table Q]] and the simulation 
results are shown in Fig. [7] to Fig. [12] 

2) 20 Users: Next we show the convergence of Algorithm 
RPD with 20 users and 64 subchannels. We also consider the 
4 cases in Table [TT] Simulation results for the convergence of 
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Fig. 9. Case 3 for 40 users: primal and dual variable convergence of 
Algorithm RPD 



Fig. 12. Case 4 for 40 users: total weighted rate convergence of Algorithm 
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Fig. 13. Case 1 for 20 users: total weighted rate convergence of Algorithm 
RPD 
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Fig. 11. Case 4 for 40 users: primal and dual variable convergence of 
Algorithm RPD 
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Fig. 14. Case 2 for 20 users: total weighted rate convergence of Algorithm 
RPD 



rates are plotted in Fig. [13] to Fig. Q~6] 
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Fig. 15. Case 3 for 20 users: total weighted rate convergence of Algorithm 
RPD 
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Fig. 16. Case 4 for 20 users: total weighted rate convergence of Algorithm 
RPD 
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3)4 Users: Next we show the convergence of Algorithm 
RPD with 4 users and 64 subchannels. Simulation results of 
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TABLE III 

Impact of System Parameters with 40 Users 



Case 
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Fig. 19. 40 users: CDF of users' achievable SNRs 

the rate convergence for the 4 cases in Table [TT] are plotted in 
Fig. [P71 and Fig. [18] Here, Cases 1-2 have the same results, 
as shown in Fig. [17] and Cases 3-4 have the same results, as 
shown in Fig. [18] 

C. Impact of System Parameters 

In Table [Till we demonstrate how different system parame- 
ters affect the convergence and performance of the algorithm. 
In all cases, we have 40 users and 64 subchannels. Users' 
channel conditions are the same as in Section [V-B II It is clear 
from Table imi that both positive self-noise and finite maximum 
SNR constraints reduce the system performance and typically 
increase the convergence time. 

In Fig. [19] we plot the Cumulative Distribution Function 
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TABLE IV 
Impact of System Parameters with 20 Users 
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TABLE V 

Impact of System Parameters with 4 Users 
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Fig. 21. 4 users: CDF of users' achievable SNRs 

(CDF) of users' allocated SNRs {^Pij^ijl (x,j + Ppijetj^j of all 
four cases after Algorithm RPD converges. In Case 2, the 
SNRs are "truncated" at 20dB due to s u = 20dB. In Case 3, the 
SNRs are also no larger than 20dB since they are effectively 
upper bounded by 1 //3. However, the effects of self-noise and 
maximum SNR constraints are apparently different. Case 4 
shows a further reduction in SNR due to both factors. 

We also give results for the case of 20 users under the same 
conditions as Section IV-B2I The simulation results are shown 
in Table HV] and Fig.gD] 

Similarly, we show results for the case of 4 users under the 
same conditions as Section IV-B3 1 in Table [V] and Fig. [2T1 



D. Comparison with the Dual-based Algorithm 

In Fig. [22] we compare the convergence speed of Algo- 
rithm RPD with the dual-based centralized optimal solution 
proposed in fl5l . The self-noise coefficient is /3 = 0.1, and 
the maximum SNR coefficient is Stj = 20dB for all i and j. 
We vary the number of users from 4 to 40. For a fixed user 
population size, we randomly generated 10 sets of different 
weights and channel conditions. We plot both the average 
and the standard deviation (i.e., error bar) of the number of 
iterations for both algorithms as the number of users changes. 
In all cases, the reduced primal-dual algorithm converges with 
less number of iterations and a much smaller variance^] 

VI. Conclusion and Future Work 

We presented the first distributed optimal primal-dual re- 
source allocation algorithm for uplink OFDM systems. The 
key features of the proposed algorithm include: (a) incor- 
porating practical OFDM system constraints such as self- 
noise and maximum SNR constraints, (b) reduced primal-dual 
algorithm which eliminates unnecessary variable updates, (c) 
distributed implementation by the end users and base station, 
(d) simple local updates with limited message passing, (e) 
global convergence despite of the existence of multiple global 
optimal solutions, and (f) fast convergence compared with the 
state-of-art centralized optimal algorithm. Currently we are 
extending this work in several directions, including proving 
the theoretical convergence speed of the algorithm. 

6 We note that the time needed for each iteration is significantly less for 
the centralized algorithm, which does not require frequent message passing 
between users and the base station. 
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Appendix 

A. Proof of Theorem [7] 

Let (x*,p*, //,*, A*) be one point satisfying the KKT con- 
dition. Motivated by the Lyapunov function used in 0, we 
consider the following La Salle function 



U v 



It is straightforward to verify that V is semi-positive definite. 
Its Lie derivative over an invariant set {(x,p,A)\ x > 0, p > 
0, A > 0} is given by 

. ^ dV ^ dV ^ dV . 



'■j 



dAi 



i.j 

+ ^ - Yj x u - ! ) + Tj {Ai - w X pv - Pi) 

j ' i j 

= ~ Sij)(Pij - P*j) 

+ Ys^j - m})( Yj x h - l) + &< - w X ph - Pi) 

j i 



•■J 



ij .. 

Pij -Pij 

X U ~ x ij 

Pij -Pij 



(29) 



where f* = fij(x^p*p and g* = gy(x*, ptX 

Now let's look at the expression in d29l ). The first sum is 
non-positive, since (A/ - gij)(pij - p*p = when py > and 
(Aj - gij)(pij - p*p < when py = according to (19[ . The 
second sum is also non-positive, since Tij x *j ^ 1 ar, d P-) — 
when Yjj x *j < 1 according to ©. In the same way, the third 
sum is also non-positive according to (TTOb . 

As to the fourth sum, it has f*. < u* according to ( [T2l . and 
the inequality holds only when x*j = 0. Similarly, g*j < A* 
according to ([T3l . and the inequality holds only when p*. = 0. 
When x*j + 0, we have p*j + and 



Hence, the fourth sum of (1291 is also non-positive. 

For the fifth sum, first we notice [fj,gij] T is the gradient 
of the following function 

PU e U 



Hijixjj, pij) = Wi{xij + eif) log 1 + 



By concavity of the function ffy, we have 
H u (x u ,pij) < Hij(x*j,Pij) + V#y(Xy,p*) 

Hij(x*j,p*p < Hij{Xij,Pij) + VHij(x u ,pij) 



Xij 


-x*. 
IJ 


Pij 


-Pij 


X ij 


— Xij 


P'j 


- Pij 



(30) 
(31) 



Consequently, the inner product of the gradient difference and 
the variable difference is non-positive, i.e. 



([fij,gij] ~ U7j,g*ji) 



■ x. . 

'j 



PU -Pij 



<0; 



hence, the fifth sum is also non-positive. 

As such, we have V < 0. According to La Salle princi- 
ple l20l . trajectories of the system in (fZTJ to (l23l converge 
to the set Vo = {(*, P, A) : V — o} globally and asymptotically. 
It is clear that at any optimal solution of Problem [2] we will 
have V = 0, i.e., such solution is in set Vo. 

Over set Vo, the five sums in ( 1291 must all equal to zero so 
as to ensure V = 0. Then we have the following observations: 

. jj.j is nonzero only if 2, x*j = 1 (according to (0 and the 
fact that the second sum in d29 i 

. Pij - p*ij 

in d29| i equals to zero). 

Combining the second observation with (l30l and (1311 , we 
further know that V(x, p,p, A) e Vo, 



([fogiA-UZ'Siji) 



equals to zero); 
= (since the fifth sum 



#y(xy,py) = Hij(x*j,p*j) + V//,,(. v ; ; , />*,-) 



X U x ij 



Pij -Pij 

Taking the derivative of both sides, we have V(x,yu, A) € Vo, 



fij(Xij,Pij) 
gij(Xij,Pij) 



VHjj(x*j,p*j) = constant. 
From gij(Xij,Pij) = gij(x*j,p*j), we can derive 



1 +Pe u - 



1 +/3e u - 



PU 



V'ij 



1 + (J3 + \) eij 



\+(j3+l)e u 



Pij 



x 'j + € U I 



1 J 



(32) 



It can be verified that the function in both sides in ( |32l , 
i.e. +/teyyj(l + (J3 + l)^,'^), is a monotonically increasing 

function of y. Therefore, Eq. d32l implies p '' = /' J . 

Suppose there exists an optimal solution {x*j,p*j} with 
Yjj p*j < Pi- For the user i on one certain subchannel /, it must 
have pa = Pi-Yjj^i P*j > P*u> which indicates that there always 
exists a new allocated value pu larger than the optimal p* By 
substituting this pu into (0 and keeping all rest p*. (Vi, Vj + I) 
unchanged, the value of the objective function in Problem|2]is 
increased. That contradicts with the assumption that {#*.,/>*.} 
is optimal. Therefore, for any optimal solution fx*., /?*.}, we 
must have 2j P*j - Pi- 

Considering the above result = (xy + eij) in set Vo, 
it has at least one nonzero py. Consequently, with the fact 
k = gij(xij, pij) according to (US and gy(xy, p u ) = g*j (which 
is always positive seen from gy's definition) discussed above, 
we get A; = A* > 0, which implies that Aj remains constant 
over set Vq. 
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B. Proof of Lemma Q] 
From (l27l i, we have 

A-iBx = P — AiBe = constant. 

Result Q] can be derived by taking derivatives (with respect to 
time) on both sides of the above equation. 

For result [2] it can be verified that the characteristic function 
of the linear system d25ll-(l26ll is a product of positive diagonal 
matrix and a skew-symmetric matrix, 







-K X A\ 




K x 
K» 





A] 







Since all eigenvalues of a skew-symmetric matrix are either 
zero or purely imaginary, the eigenvalues of the characteristic 
function are also imaginary. Hence trajectories of the linear 
system do not converge. 

C. Proof of Theorem [2] 

By linear system theory, fi is completely observable from 
the constant AiBx if and only if the complete observability 
(|28T > holds ETI . Then the invariant set Vo contains only the 
equilibria of the linear system in d25l l to (l27l i. which are 
the global optimal solutions of Problem [2] Consequently, all 
trajectories of Algorithm RPD converge globally and asymp- 
totically to the global optimal solutions. 



D. Proof of Corollary Q] 

Let us choose matrix K x to be kl where k is a positive 
constant. By direct computation, we get A2BK x A T l = kA2BA^. 

First consider the case where N < M, it can be observed 
that A2BK X A\ (the upper part of the matrix in d28l l) has rank 
N, because A2BA\ has N linearly independent columns under 
such condition. 

Now consider the case where M < N, in which case the rank 
of matrix kAiBA\ will be less than N. In this case we need to 
look at matrix kK^AiA^ — crl (the bottom part of the matrix in 
(ED). Consider the diagonal matrix K> l A l K x Al = kK»A x A\, 
which can be written as 
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kKf 



M 



If all diagonal elements of matrix take different values, the 
rank of matrix kK^A\A T y - crl would be N - 1, because the 
eigenvalue cr is one of the diagonal elements of the diagonal 
matrix kK^AiA^. Combining with A 2 BK X A^ which has full 
row rank M > 1, the combined matrix in (f28b satisfies the 
observability condition in Theorem [2] To sum up, when K x = 
kl and all update stepsizes kfj (j = 1 , . . . , N) take different 
values, the observability of the linear system is guaranteed. 



